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ABSTRACT 


A  locally-implicit  scheme  for  steady-state  solution  of  the  thin-layer  Navier- 
Stokes  Equations  is  simplified  by  elimination  of  coefficient  matrices.  This  reduces 
both  arithmetic  computation  and  computer  storage  requirements.  An  added 
benefit  is  the  simplification  of  the  algorithm  which  eases  the  coding  task.  The 
locally-implicit  scheme  uses  finite-volume  spatial  discretization,  locally-implicit 
time  integration,  Jameson-type  artificial  dissipation  terms  and  a  modified  Gauss- 
Seidel  iteration.  The  modified  method  is  tested  for  subsonic  and  transonic  flows 
over  an  RAE  2822  airfoil. 


111 


TABLE  OF  CONTENTS 


CHAPTER  PAGE 

1.  INTRODUCTION  . 1 

2.  LOCALLY-IMPLICIT  SCHEME  FOR 

ONE-DIMENSIONAL  MODEL  EQUATION  . 3 

2.1  Model  Equations . 3 

2.2  Euler  Implicit  Scheme  for  Model  Equation . 4 

2.3  Local  Stability  Analysis  for  Model  Equation . 8 

3.  LOCALLY-IMPLICIT  SCHEME  FOR 

NAVIER-STOKES  EQUATIONS  . 14 

3.1  Original  Locally-Implicit  Scheme  for  Navier-Stokes  Equations  .  14 

3.2  Modification  of  Locally-Implicit  Scheme  for 

Navier-Stokes  Equations  . 19 

4.  RESULTS . 21 

5.  CONCLUDING  REMARKS . 30 

LIST  OF  REFERENCES . 31 

APPENDIXES . 33 

APPENDIX  1.  DERIVATION  OF  Cjt k . 34 

APPENDIX  2.  ADDITIONAL  SMOOTHING  FACTORS 

IN  SCHEME . 41 

VITA . 45 


IV 


LIST  OF  FIGURES 


FIGURE  PAGE 

1.  Stability  Plot  of  Modified  Locally-Implicit  Scheme  Applied  to 

1-D  Equation  (3),  (u  =  . 13 

2.  Grid  for  RAE  2822  Airfoil  (128  x  32  Cells)  . 23  ' 


3.  Results  for  Case  2  (Moo  =  0.676,  a  =  2.40  deg.,  Re  =  5.7xl06)  .  .  24 

4.  Results  for  Case  4  (Moo  =  0.725,  a  =  2.92  deg.,  Re  =  6.5x10s)  .  .  26 

5.  Results  for  Case  5  (Moo  =  0.730,  a  —  3.19  deg.,  Re  =  2.7xl06)  .  .  28 


LIST  OF  SYMBOLS  AND  ABBREVIATIONS 


a  see  Equation  (2),  also  speed  of  sound 

A,  B  Jacobian  matrices  corresponding  to  inviscid  flux  vectors 

c  see  Equation  (6) 

CFL  Courant  number 

Cp  pressure  coefficient 

d  artificial  dissipation 

du  correction  quantity  for  A Uj  [see  Equation(8)] 

e  energy 

E ,  F  inviscid  flux  vectors 

g  amplification  factor 

i  \/~T 

I  identity  matrix 

J  Jacobian  of  metric  tensor  for  coordinate  transformation 

L  any  operator 

M  Mach  number 

p  pressure 

Pr  Prandtl  number 

Q  vector  of  dependent,  conservation  variables 

A Q,  A u  change  in  the  solution  per  time  step 

r  see  Equation(6) 

R  simplified  viscous  flux  Jacobian  (diagonal)  matrix 

Re  Reynolds  number 

Res  residual 

RHS  right  hand  side 


VI 


s 

viscous  flux  vector 

t 

time 

u 

general  one-dimensional  variable,  Cartesian  velocity  component 

in  x  direction 

V 

Cartesian  velocity  component  in  y  direction 

a 

see  Equation  (20),  also  angle  of  attack  (see  Table  1) 

7 

ratio  of  specific  heats 

e(2) 

factor  used  in  artificial  dissipation 

£<4> 

factor  used  in  artificial  dissipation 

C 

see  Equation  (20) 

K<2) 

second- difference  dissipation  constant 

k<4) 

fourth-difference  dissipation  constant 

P 

dynamic  viscosity 

V 

factor  used  in  artificial  dissipation,  also  kinematic  viscosity 

p 

density 

T 

time 

U 

relaxation  factor 

Subscrip 

ts 

j 

x- direction 

k 

y-direction 

til 1 

curvilinear  coordinates 

oo 

free  stream  condition 

Superscr 

ipts 

m 

iteration  sweep  count 

n 

time  level 

vii 

CHAPTER  1 


INTRODUCTION 

Though  computational  fluid  dynamics  (CFD)  has  been  in  existence  for  ’ess 
than  three  decades,  it  is  already  recognized  as  an  increasingly  powerful  tool  for 
aerodynamic  design  of  aerospace  vehicles1.  Flight  test,  ground  test,  and  CFD 
compliment  one  another  as  aerospace  vehicle  design  requirements  become  more 
demanding.  Besides  being  used  directly  in  the  design  process,  CFD  also  plays  an 
important  role  in  furthering  scientific  understanding  of  complex  flow  phenomena. 

Within  the  CFD  spectrum  are  various  levels  of  maturity.  At  one  end  of  the 
spectrum  are  panel  codes  that  solve  linear-type  flows  for  which  CFD  is  quite  ma¬ 
ture.  These  codes  can  readily  handle  very  complicated  geometries  in  reasonable 
computer  run  times.  At  the  other  end  of  the  spectrum  are  the  Euler  and  Navier- 
Stokes  codes  which  handle  the  complex  physics  of  nonlinear-type  flow.  These 
codes  currently  require  substantial  computer  resources  (both  memory  and  time) 
and  are  limited  to  simple  geometries  (relative  to  linear- type  flow  solvers),  but 
provide  essential  physical  insight  into  many  classes  of  flow  problems.  To  make 
Euler  and  Navier-Stokes  codes  even  more  useful  than  they  currently  are,  gains 
must  be  made  in  several  areas  such  as  improvement  in  speed  of  convergence; 
improvement  in  handling  complex  geometry;  reduction  in  numerical  diffusion 
for  flows  containing  shocks,  wakes  and  vortex  structures;  and  better  turbulence 
modeling1.  The  scheme  discussed  in  this  thesis  addresses  the  area  of  increasing 
computational  efficiency  of  algorithms. 

Current  useful  Euler  and  Navier-Stokes  schemes  fail  into  one  of  two  cate¬ 
gories:  explicit  or  implicit  methods.  MacCormack2,  in  comparing  these  meth¬ 
ods,  points  out  the  need  to  use  an  implicit  procedure  to  avoid  restricting  CFL 


numbers  to  small  values.  He  states  that  Newton  iteration  and  line  Gauss- 
Seidel  procedures  show  potential  for  significantly  increasing  numerical  efficiency. 
Chakravarthy3  has  shown  various  relaxation  procedures  for  implicit  schemes. 
Reddy  and  Jacocks4  have  prerented  a  locally-implicit  scheme  for  solving  the 
Euler  equations  for  large  problems  using  a  modified  Gauss-Seidel  relaxation  pro¬ 
cedure.  They  find  a  CFL  number  of  about  ten  to  be  appropriate  for  their  scheme. 
Reddy  and  Nayani5,  and  Nayani6  have  extended  the  locally-implicit  scheme  to 
solution  of  the  two-dimensional  Navier-Stokes  equations.  The  work  presented  in 
this  thesis  is  a  modification  to  that  scheme.  The  modification  results  in  elimina¬ 
tion  of  seven  four-by-four  matrices  thereby  making  the  scheme  matrix-free.  The 
modified  scheme  involves  less  arithmetic  computation  and  requires  less  computer 
storage.  Primary  motivation  for  the  work  presented  in  this  thesis  is  to  establish 
the  viability  of  the  matrix-free  scheme. 

The  scheme  Nayani  presented  in  his  dissertation6  is  referred  to  throughout 
this  thesis  as  the  “original”  scheme.  The  modification  reported  in  this  thesis 
is  referred  to  as  the  “modified”  scheme.  Nayani ’s  computer  code  was  used  as 
the  starting  point  for  the  results  presented  in  this  thesis.  Parametric  studies 
of  various  relaxation  parameters  with  regard  to  the  convergence  process  and 
multigrid  acceleration  of  convergence  of  the  scheme  are  not  addressed  here,  since 
such  features  would  be  similar  to  the  original  scheme. 
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CHAPTER  2 


LOCALLY-IMPLICIT  SCHEME  FOR  ONE-DIMENSIONAL 
MODEL  EQUATION 

In  this  chapter,  a  model  equation  is  used  to  describe  the  locally-implicit 
method.  It  starts  with  a  finite  difference  scheme  based  on  central  differences  for 
spatial  derivatives  and  Euler  implicit  time  integration.  A  modified  symmetric 
Gauss-Seidel  iteration  is  developed  for  solving  the  implicit  equations  at  every 
time  step.  The  scheme  is  shown  to  be  unconditionally  stable  by  Fourier  stability 
analysis. 


2.1  Model  Equations 

Nayani6  has  demonstrated  the  basic  locally-implicit  scheme  and  its  stability 
with  the  one-dimensional  diffusion  equation: 


du  _  d2u 
~dt~ud^ 


(1) 


Reddy  and  Nayani5  have  established  the  unconditional  stability  of  the  locally- 
implicit  scheme  for  the  linearized  Burgers  equation, 


du  du  _  d2u 
dt  +  °  dx  V  dx2 


(2) 


which  serves  as  a  model  for  the  Navier-Stokes  equations.  The  coefficient  a  rep¬ 
resents  a  convection  velocity  which  is  independent  of  x  and  can  be  positive  or 
negative.  The  coefficient  v  represents  a  generalized  diffusion  coefficient  corre¬ 
sponding  to  the  fluid  viscosity  or  an  artificial  dissipation  coefficient  designed 
to  be  effective  near  shocks  ( v  «(2)|a|A:r,  where  is  a  second- difference 
dissipation  constant  obtained  through  numerical  experimentation). 
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The  following  one-dimensional  model  equation  will  be  used  in  this  thesis  to 
present  the  modified  locally-implicit  scheme  and  its  stability. 


du  du  d*u 

dt  "I"  Q  dx  U  dx4 


(3) 


In  Equation  (3),  v  represents  an  artificial  dissipation  coefficient  designed  to 
be  effective  everywhere  in  the  flow,  except  near  shocks  (v  %  K(4,|a|  Ax3,  where 
is  a  fourth-difference  dissipation  constant).  Artificial  dissipation  of  this  type 
has  been  shown  by  Jameson7  and  others  to  suppress  the  non-linear  instabilities 
which  arise  in  central  difference  schemes  for  convection  dominated  flows. 


2.2  Euler  Implicit  Scheme  for  Model  Equation 

Reddy  and  Nayani5  have  analyzed  the  locally-implicit  scheme  for  Equation 
(2)  with  central  difference  approximations  for  the  spatial  derivatives  and  Euler 
implicit  time  integration. 


..n+i 

i  ui- 1  _  ui+ 1 
2Ax 


-  2u"+1 


+  “  K 


The  same  basic  scheme  will  be  demonstrated  for  Equation  (3).  As  in  Equation 
(4),  the  Euler  implicit  time  integration  and  central  spatial  difference  approxima¬ 
tions  for  the  spatial  derivatives  are  used. 


-  u"4,1  un+i 

■7-1  _  v  J  +  2 


-  4 +  6uy +1  -  4 +  un+l 


A  direct  solution  of  Equation  (5)  along  with  its  accompanying  boundary  con¬ 
ditions  requires  solution  of  a  pentadiagonal  system  of  equations  for  each  time 
step.  In  multi  dimensions,  the  matrix  is  too  big  for  a  direct-solution  method. 
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Therefore,  the  locally-implicit  method  described  by  Reddy  and  Nayani5  is  used 
to  obtain  an  asymptotic  steady-state  solution  of  Equation  (3). 

The  delta  form  of  Equation  (5)  is  sought  by  rewriting  Equation  (5)  with  the 
use  of  the  following  definitions  for  Aw,  c  and  r. 


Auj  =  u™+1  —  u?,  c  =  =  CFL,  r  = 

J  J  J  At 


aAt 


_  vAt  «^4^|a|At 


Ax 


Ax4 


Ax 


=  /c^lcl 


c 

Au_,-  +  1  -  Auj-i) 

+  r(Auj+2  —  4Auj+i  +  6Auj  —  4Au_?_1  +  Auj-2 )  =  Res™  (6) 

Res™  =  -|(«7+i  -  -  r(u™+2  -  4 u]+1  +  6u™  -  4 u?_2  +  u™_2 )  (7) 

Res”  is  the  residual.  A  steady-state  solution  is  obtained  for  Equation  (5)  by 
driving  the  residual  to  zero. 

In  the  standard  Gauss-Seidel  iteration  procedure,  when  sweeping  left  to 
right  ( j  increasing),  Auj  is  solved  for  by  setting  AuJ+j  to  zero  and  using  the  value 
for  Auj-i  available  from  the  last  iteration  step.  This  procedure,  however,  is  not 
stable  because  of  the  convection  terms,  even  if  a  symmetric  Gauss-Seidel  scheme 
is  used.  Therefore,  a  modified  Gauss-Seidel  procedure4  that  is  unconditionally 
stable  is  used  which  employs  an  inner  iteration.  First,  denote 

duj  =  Aujm+1^  —  A 


where  m  is  the  inner  iteration  index  count  and  duj  is  the  correction  quantity  for 
Aiij  in  the  mth  inner  iteration  sweep.  Substituting  this  expression  for  duj  into 
Equation  (6)  gives 


(8) 


*  c 

duj  +  ~(duj+ 1  -  duj-i) 

+  r(duj+2  —  4  duj+i  +  6  duj  —  iduj-i  4-  duj- 2)  =  -RHS 

Bi/S  =  Res]  -  |Au<m)  +  |(Aix£J  -  Au’">) 

+  r(AU<™>  -  4Au<">  +  6Auf>  -  4Au<™>  +  Au$”>)] 

For  a  left-to-right  sweep,  the  newly  computed  du's  (i.e.,  duj- 1  and  du.j-2)  are 
brought  over  to  the  right-hand  side  of  Equation  (8)  producing  the  following 
equation. 


duj  -f  ~(duj+1)  +  r(duj+2  —  4duj+i  4-  6duj)  =  RHS 


(9) 


RHS  =  Res]  -  Lj(  Au) 
ly(Au)  =  AU<m|  +  §(Au£>  -  Au^+])) 

+  r(Au<  -  4Au£>  +  6Auf>  -  4Au£+1)  +  Au^1*) 

So  fax,  this  development  is  identical  to  that  shown  by  Nayani6  for  the  heat 
equation  and  is  the  departure  point  for  the  present  scheme. 

The  latest  available  values  for  uj  are  denoted  by  starred  quantities  as  follows. 


=  u? 


*  (m) 


Uj_2  and  Uj+ 2  have  similar  representations.  Similarly,  in  the  right-to-left  sweep, 
starred  quantities  represent  the  latest  available  values  also.  The  terms  in  the 
residual  [Equation  (7)]  can  be  combined  with  all  the  other  terms  in  RHS  except 
A u^m\  to  produce 


(11) 


rhs  =  -  A«$m) 

-  r(u*+2  -  4u}+1  +  6 u*j  -  4 u*_^  +  u*_2 ) 


or 

(12) 
(13) 


duj  +  -(<2ttj+i)  +  r(duj+ 2  —  4duJ+i  +  6duj)  =  —  Au^m)  (14) 

If  the  right-hand  side  is  driven  to  zero,  a  time-accurate  solution  to  the  discrete 
unsteady  problem  [Equation  (5)]  for  u"+1  results.  The  left  side  of  Equation  (14) 
can  be  modified  to  achieve  efficiency  and  convergence  of  the  scheme  without 
altering  the  actual  solution. 

First,  approximate  duj+2  and  duj+i  by  duj.  This  gives 

(1  +  ^  +  3 r)duj  =  Res*  -  A (15) 
Using  this  approximation,  for  a  right-to-left  sweep,  Equation  (14)  is 

(1  —  ^4-  3 r)duj  =  Res*  -  A Uym)  (16) 

This  type  of  symmetric  Gauss-Seidel  iteration  is  not  stable.  The  symmetric 
Gauss-Seidel  iteration  is  modified  by  replacing  c  by  |c|  and  using 


RHS  =  Res*  -  A u<m) 

Resj  =  -^(ui+i  “  uj- 1)  “  r(uj+ 2  ~  *uj+  i  +  -  4u*_x  +  u*_2) 

where  ReSj  is  the  residual  computed  for  the  latest  values  of  u. 

The  left-to-right  scheme  is  given  as 
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(1  +  Y  +  3  r) 


as  the  coefficient  of  duj  for  both  left-to-right  and  right-to-left  sweeps. 

For  initialization  of  the  inner  iteration,  set  Au|0)  =  0.  Thus,  the  modified 
scheme  is  represented  by 

(1  +  ^  +  3r)duj  =  Res*  -  Au<m)  (17) 

None  of  the  modifications  have  jeopardized  the  ability  to  obtain  time- 
accurate  solutions  to  Equation  (5).  For  time  accuracy,  symmetric  sweeps  per 
time  step  are  performed  until  the  right-hand  side  of  Equation  (17)  becomes  suf¬ 
ficiently  close  to  zero.  A  symmetric  sweep  for  this  one- dimensional  case  is  one 
left-to-right  sweep  followed  by  one  right-to-left  sweep.  In  this  thesis,  since  the 
steady-state  solution  is  sought,  only  one  or  two  symmetric  sweeps  per  time  step 
axe  necessary  while  advancing  the  solution  in  time  until  desired  convergence  is 
obtained  (i.e.,  until  Res*-  is  driven  to  zero). 

The  major  advantage  of  this  modification  to  the  scheme  is  the  resulting 
reduction  in  both,  arithmetic  computation  and  computer  storage  requirements. 
This  fact  will  become  evident  in  Section  3.2. 

2.3  Local  Stability  Analysis  for  Model  Equation 

Local  linear  stability  analyses  for  this  locally-implicit  scheme  have  already 
been  presented  for  two  of  the  three  model  equations  of  interest.  Both,  Burgers’ 
equation5  and  the  heat  diffusion  equation6  have  been  shown  to  be  uncondition¬ 
ally  stable  in  a  local  linear  sense  for  all  CFL  numbers.  The  stability  analysis  for 
the  third  model  equation  of  interest  [Equation  (3)]  is  presented  in  this  section. 
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The  standard  Fourier  stability  analysis  is  performed  on  Equation  (17)  after 


expanding  Res*j  from  its  *  notation  back  to  the  pure  delta  form. 


(1  +  ^  +  Zr)duj  =  RHS  (18) 

RHS  =  RHSL-r  =  Res J  -  (Au'm)  +  |(Au$?>  -  Au$"+1)) 

+  r(AuJ”J  -  4AU<">  +  6Au'm)  -  4Au'”+1)  +  Au^2+1>)] 


or 


RHS  =  RHSr-L  =  Res J  -  [Au'”'  +  |(At4+,+1)  -  Au$™|) 

+  -  4Au£1+,>  +  6Au<m)  -  4AuJ”>  +  Au<">)] 

Res "  =  -[|(“”+i  -  «"-i)  +  K“"+2  -  4u”+1  +  6u”  -  4u"„,  +  u"_2)] 

du,  =  Au'"+1)  -  Au$m| 

One  symmetric  sweep  will  be  considered  in  this  analysis.  For  the  left-to- 
right  sweep,  m  =  0.  Recalling  that  Au^  is  set  to  zero,  Equation  (18)  becomes 

(1  +  "  +  3  r)duj  =  RHSL-R  (19) 

RHSl-r  =  Res] J  -  [|(-Au$1_,1)  +  r(-4Au;i_)1  +  Au^)] 

dUj  =  AuJ1}  -  Au$0)  =  Au^ 
uj  =  u”  +  Au$1} 
u]+1  =  unj  +  A  U(P 
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Discrete  modal  solution  of  Equation  (18)  is  sought  in  the  form 


un  =  VneijQAx  =  Vnetj<,  i  =  -/-l,  C  =  <*Ax,  0  <  C  <  tt 

un  =  yheij  ( 

A u<3)  =  uj  -  u]  =  ( Vh  -  Vn)eiK 
A uf  =  u]+1  -  unj  =  {Vn+1  -  Vn)eijc 

Substituting  Equations  (20)-(22)  into  Equation  (19)  and  dividing  each 
the  resulting  equation  by  e,J<*  and  then  simplifying  gives 


Ti(V"  -  Vn)  =  T2Vn 
b  =  1  +  —  +  3r 

2\  =  6- [|  +  r(4-e "*)]«"* 

J2  =  -icSinC  -  r[2Cos(2()  -  8 Cos(  +  6] 

The  amplification  factor  g  for  the  scheme  is  defined  by  the  equation 


Vn+1  =  gVn 


So,  for  the  left-to-right  sweep  the  amplification  factor  g  is 


-  V *  ,  T2 
9  Vn  ~  1  +  Tx 

For  the  right-to-left  sweep,  Equation  (18)  becomes 


(l  +  y  +  3  r)duj  =  RHSR.L 


(20) 

(21) 

(22) 

(23) 

side  of 

(24) 


(25) 


(26) 


(27) 
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RHSR-L  =  Res ?  -  [A +  |(A -  At#?,) 

+  r(Ait^2  -  4Au^h\  +  6Au^1}  -  4Au^1_)1  +  Au^)] 
duj  —  Auy2^  —  Au^ 

Equations  (20)-(23)  are  substituted  into  Equation  (27)  and  each  side  of  the 
resulting  equation  is  divided  by  e'-?l*  and  then  simplified  giving 


T3(vn+1  -  vn)  =  T2vn  +  r4(vft  -  Vn) 


(28) 


T3  =  b  +  ll  -  r(4  -  e‘<)]e^ 

T,  =  M  -  3r  +  [|  +  r(4  -  e -*))«-« 

After  the  right-to-left  sweep  is  performed  the  amplification  factor  for  the  scheme 


is 


yn+l  rp  T 

f=  — =  1  +  5f  +  5f(j-l)  (29) 

Substituting  Equation  (26)  into  Equation  (29)  and  simplifying  gives  the  scheme 
amplification  factor  as 


9  =  1  +  jr(l  +  jr)  (30) 

The  stability  requirement  is 

l?l  <  1  (31) 

Therefore,  to  have  a  stable  scheme,  the  following  must  be  true. 

b|  =  |l  + J(l  +  ^)|<1  (32) 
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Figure  1  shows  a  plot  of  the  modulus  of  the  amplification  factor  versus  the 
phase  angle  for  Equation  (32).  CFL  varies  from  one  to  100  and  u  =  As 
stated  at  the  beginning  of  this  section,  the  locally-implicit  scheme  is  indeed  stable 
for  the  one-dimensional  model  equation  which  models  the  fourth-order  numerical 
dissipation.  From  the  other  stability  analyses  that  have  been  presented  by  Reddy 
and  Nayani5,  and  Nayani6,  a  CFL  number  of  10  was  chosen.  This  is  a  good  choice 
for  this  model  equation  as  well. 


2 


Figure  1.  Stability  Plot  of  Modified  Locally-Implicit  Scheme  Applied 
to  1-D  Equation  (3),  (u  = 
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CHAPTER  3 


LOCALLY-IMPLICIT  SCHEME  FOR  NAVIER-STOKES  EQUATIONS 

Nayani  gives  a  detailed  derivation  of  the  original  scheme  in  his  dissertation6. 
Relevant  equations  of  the  scheme  are  presented  in  Section  3.1  without  rederiva¬ 
tion.  The  technique  he  has  presented  possesses  the  following  salient  features: 
finite  volume  spatial  descretization,  Jameson- type  artificial  dissipation  terms7, 
modified  Gauss-Seidel  inner  iteration.  For  the  steady-state  solution,  he  has  used 
local  time  integration  which  varies  the  time  step  cell-by-cell  in  order  to  maintain 
a  constant  CFL  number  throughout  the  computational  grid.  Multigrid  technique 
has  been  used  for  acceleration  of  convergence,  but  that  aspect  of  the  code  is  not 
addressed  in  this  thesis.  He  has  applied  the  scheme  to  the  Navier- Stokes  Equa¬ 
tions  incorporating  the  thin-layer  approximation  and  using  the  Baldwin- Lomax8 
turbulence  model.  The  modification  to  the  original  scheme  (presented  in  section 
3.2)  does  not  change  any  of  the  above  features. 

3.1  Original  Locally- Implicit  Scheme  for  Navier-Stokes  Equations 

Nayani6  starts  his  development  with  the  Navier-Stokes  Equations  in  nondi- 
mensional  form  for  generalized  curvilinear  coordinates.  Then  he  makes  the  thin- 
layer  approximation  and  simplifies  the  equations  to  arrive  at 


+  0-(-y(E  +  Z(F)-R'-'±S=  0 


(33) 
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where 


'  p  ' 

pu 

pu 

pu 

,  E  = 

pu 2  +  p 

,  F  = 

puv 

pv 

puv 

pu  +  P 

.  e  . 

.  u(e  +  p)  . 

,v(e  +  p)  . 

S=J 


-l 


{Vxmi  +  Vym2) 

.  (r)xm2  +  riym3) 

L  ■qx(um1  +  vm2  +  m4)  +  rjy(um2  +  vm3  +  m5 ) 


mi  =  |(4 VxUv  -  2 TjyVr,) 
m2  =  n(riyUr,  +  r)xvv ) 
m3  =  ^(-2t)xuv  +  4r)yvv) 

m4  =  Pr(y-  1)^(q2) 
ms  =  PrT7~i')T?va,?(a2) 

After  integrating  the  Navier-Stokes  Equations  using  a  finite-volume  ap¬ 
proach,  he  derives  the  Euler  implicit  scheme  in  delta  form  using  the  modified 
Gauss-Seidel  iteration  technique. 


<34) 

wout  is  an  outer  iteration  relaxation  parameter.  AQ  is  produced  from  the  inner 
iteration 


A  Q^k+1)  =  A  gj.y  +  windg^k,  m  =  0, 1, 2, 3  (35) 

u>in  is  an  inner  iteration  relaxation  parameter,  m  is  the  sweep  number.  Four 
sweeps  (identified  as  m  =  0, 1,2,3)  using  Equation  (35)  make  up  one  symmetric 


sweep  for  this  two-dimensional  case.  The  A Q^l  term  is  then  used  as  A Qjk  in 
Equation  (34)  which  is  added  to  Q™ k  to  give  the  new  value  Q^1  ■ 

The  values  for  dQJtk  are  produced  from 


Cjtk  ■  dQjj  =  Resj  k  -  Ljik(AQ) 


Resl„  =  -[!/,£-  -  -  [-«£”  + 

h~t<k 


jyk-% 


+  Re~\Sn) 


j,k+i 


— c?  I  +  d 
j )  k  —  j1  j  ^ » k 


j,k+± 


j'k-i 


(36) 


(37) 


Re  is  the  Reynolds  number,  d  is  the  artificial  dissipation  flux  of  Jameson7  type 
given  as  follows. 


~e\4+l  k(Qj+2,k  ~  3^>+l,*  +  3 Qj,k  ~  Qj-l,k) 


>  +  *,* 


where 


=  ^(2)  ’  vhk) 

^4+)i,fc  =  Max[0,(K(4)-eJ2+)iifc)] 


vi,k  = 


lPi+i.fc  ~  2Pi,k  +Pj-i,fcl 
|p;+i,fc  +  2Pj,k+Pj-i,k 
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d  is  defined  similarly  for  other  indices.  C  hk  is  a  diagonal  matrix  which  is  defined 
from  a  heuristic  derivation  (see  Appendix  1). 


=  CFI(g)(I  +  ^)/+ 

+ +  %*•**.»' + 3%^- id 


At 

r— 1 


+  3^±!±e(4)  j 

At  J>fc+s  At  J>fc+2 


r-i 

At  V'i.fc-T  ‘  "  At 


(38) 


where  i?  is  a  diagonal  matrix  derived  from  the  approximation  of  (see  Ref¬ 
erence  6).  Note  that  Nayani’s  dissertation6  contains  a  misprint  for  Cy,*  [his 
Equation  (68)].  Equation  (38)  above  reflects  the  correct  Cjyk-  The  operator  Lhk 
of  Equation  (36)  is  similar  to  Lj  found  in  Equation  (9)  for  the  one-dimensional 
case. 


+  JJg±e%hkAQi+2,t 

J-!  J- 1 

I  !  ,(4)  A  |  (4)  A 

+  ~A ”  AQ'’*"2  +  ^AQj’k+2 

+  CLjtk&Qj-\,k  +  CBj'k&Qj,k- 1  +  CCj^AQj,k 

+  CRjtkAQj+i,k  +  CTjtkAQjtk+i 


(39) 


where  AQ’s  are  the  latest  available  values  (i.e.,  some  AQ’s  are  at  the  sweep 
level  (m)  and  some  are  at  the  level  (m  +  1),  depending  on  the  particular  sweep 
direction  at  the  time  Equation  (39)  is  being  applied).  The  coefficient  matrices 
CCjtk,  CLjtk,  CRj,k,  CBjtk,  and  CTjtk  are 
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CCj.fc  —  ^[(yvAn  —  xvBn)j+±j.  —  ( yvAn  —  x, 1Bn)j_^k] 

+  2^~ytAn  +  xtBn)j,k+$  -  (-y«^n  + 


+  Re  \Rjk+L+R.k_±)  +  -J±l 


,  Jj±±(2)  j  2Jj+>’k  r(4)  /4-iziiiJ2)  T 

At  i+2'fc  At  i+2>*  At  ei~bk 

J-1  7-1  7-1 

,  3  M’*  (4)  r  ,  i-*+i  (2)  j  „  M+i  (4)  j 

+  3  Ar  V},*J  +  — VHf+3~V+j; 

+  %if',l1;  +  S-T:ii14! 

At  2  At  2 


T  —  1 

—  _i(S,  4n  _  X  .  _  ->+2'fc  (4)  r 

T“ 1  J  — 1 

M'fc  (2)  J  o  J~M.  (4)  r 

At  >~2’*  At  /-i>* 


V 

J"1  J”1 

_  o  J~*~2’*  .(4)  7-  J  .(4)  r 

At  J+2’k  At  i~hk 


CBi,t  =  -^-V(A"  +  -  Re-'R^-l 


J-1 

J,*+i 

e(4)  J  - 

J-.1  , 

-(2)  f 

^7*.i 

0  2 

At 

At  >-fc-2 

3  Ar 

CTj,k  —  ( 


M±j  r<2)  T  o  J»fc+j.(4)  7-  J.*-*  (4)  r 

At  >>*+2  At  -7'*+ 2  At  hk~V 
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where 


A 


dE  dF 

dQ'  ~  dQ 


The  above  original  scheme  is  designed  to  compute  the  steady-state  solution 
to  the  thin-layer  Navier-Stokes  Equations  (33).  The  solution  is  obtained  when 
Res^k  [Equation  (37)]  is  near  enough  to  zero. 

The  Ljtk  operator  reported  by  Nayani6  is  given  by  Equation  (39).  His  code 
contains  some  additional  smoothing  factors  that  do  not  appear  in  his  disser¬ 
tation.  The  smoothing  factors  make  the  code  more  robust.  For  details,  see 
Appendix  2. 


3.2  Modification  of  Locallv-Implicit  Scheme  for  Navier-Stokes  Equations 

The  modification  to  the  original  scheme  eliminates  the  need  for  the  coef¬ 
ficient  matrices  of  Equations  (40)-(44).  By  applying  the  procedure  described 
in  Section  2.2  for  the  one-dimensional  equation  to  the  two-dimensional  Navier- 
Stokes  equations,  the  resulting  scheme  becomes 


C j,k  •  dQjtk  —  ReSj  k  —  AQ ; k 


At 


(45) 


AC?' t  =  Q‘,k  -  Ql„ 


(46) 


= -[»,£•  -  *,f”ir+!'‘  -  [-»{£* + *<*••] 

ij-f,* 


J.k+j 

j,k-j 


+  Re~1(S*) 


j+$,k 

~f~  d 

,  +d * 

j,k-k 

j~\ ,k 

J,k+  2 

;4- 2 


(47) 
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All  but  one  of  the  terms  in  LJtfc(AQ)  [Equation  (39)]  combine  with  the  terms  of 
Res™  k  [Equation  (37)]  to  produce  Res *  fc,  which,  analogous  to  the  one-dimensional 
scheme,  is  the  residual  computed  from  the  most  recent  values  of  Q  [some  values 
are  at  the  sweep  level  (m)  and  some  are  at  the  sweep  level  (m  +  1)]. 

The  main  advantage  of  this  formulation  of  the  scheme  is  the  elimination  of 
the  four-by-four  coefficient  matrices  CCj, jt,  CLjtk ,  CRjt *,  CBjtk ,  CTj [Equa¬ 
tions  (40)-(44)]  and  the  four-by-four  flux  Jacobian  matrices  An  and  Bn  (note 
that  An  and  Bn  require  significant  computer  memory  since  they  are  stored  at 
every  node).  This  eliminates  all  of  the  four-by-four  matrices  thereby  making 
the  scheme  matrix  free.  These  would  be  five-by-five  matrices  for  the  three- 
dimensional  Navier-Stokes  equations.  For  a  two-dimensional  flow  problem  of 
128  x  32  cells,  this  modification  reduces  arithmetic  computation  and  memory 
storage  requirements  by  approximately  30  percent  and  50  percent,  respectively. 
The  task  of  coding  the  scheme  is  simplified  as  well. 
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CHAPTER  4 


RESULTS 

Three  of  the  test  cases  Nayani  used  have  been  computed  using  the  modified 
scheme  and  are  compared  with  the  original  scheme  and  experimental  data9.  The 
cases  are  shown  in  Table  1  below  for  flow  over  the  RAE  2822  airfoil. 


Table  1.  Test  Cases  for  Modified  Scheme 


CASE 

M0 o 

a(deg) 

Re 

2 

0.676 

2.40 

5.7xl06 

4 

0.725 

2.92 

6.5xl06 

5 

0.730 

3.19 

2.7xl06 

The  grid  for  the  RAE  2822  airfoil  is  shown  in  Figure  2.  It  is  exactly  the 
same  grid  Nayani6  used.  An  algebraic  turbulence  model  (eddy  viscosity  model 
developed  by  Baldwin  and  Lomax8)  has  been  used  to  obtain  effective  dynamic 
viscosity  values  at  each  cell.  The  results  are  shown  in  Figures  3-5. 

Static  pressure  contours  for  the  three  cases  using  the  modified  scheme  are 
shown  in  Figures  3a,  4a,  and  5a.  The  pressure  coefficient  plots  (Figures  3b, 
4b,  and  5b)  show  that  the  original  and  modified  scheme  give  the  same  results. 
Case-by-case  comparisons  of  convergence  behavior  for  the  original  and  modified 
schemes  are  shown  in  Figures  3c,  4c,  and  5c.  Those  too  are  practically  identical. 

For  all  the  computations,  the  CFL  number  has  been  varied  over  the  first 
250  iterations  since  impulsive-start  boundary  conditions  on  the  airfoil  have  been 


21 


used.  For  the  first  150  iterations,  the  CFL  number  is  set  to  one.  For  the  next 
100  iterations,  it  is  set  to  five.  For  the  remaining  iterations,  the  CFL  number 
is  10.  There  are  alternate  ways  of  starting  the  solution  process,  such  as  gradual 
implementation  of  boundary  conditions  over  a  number  of  time  steps  instead  of 
impulsive-start,  but  such  techniques  are  not  demonstrated  here. 
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Figure  2.  Grid  for  RAE  2822  Airfoil  (128  x  32  Cells) 
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O  Experiment 
A  Numerical  (orig) 
o  Numerical  (mod) 


Figure  3.  (Continued) 
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b.  Pressure  Coefficient 


c.  Convergence  History 


Figure  4.  (Continued) 
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O  Experiment 
A  Numerical  (orif) 
o  Numerical  (mod) 


x/c 


b.  Pressure  Coefficient 


c.  Convergence  History 


Figure  5.  (Continued) 
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CHAPTER  5 


CONCLUDING  REMARKS 

The  locally-implicit  scheme  developed  by  Reddy  and  Jacocks4  which  Reddy 
and  Nayani5,  and  Nayani6  presented  for  the  two  dimensional  Navier-Stokes  equa¬ 
tions  has  been  modified  thereby  eliminating  the  need  to  compute  seven  four-by- 
four  coefficient  matrices  ,  of  which  two  axe  stored  at  every  node  of  the  computa¬ 
tional  mesh.  This  more  elegant  matrix-free  representation  of  the  locally-implicit 
scheme  results  in  the  reduction  of  both,  arithmetic  computation  and  computer 
memory  storage  requirements.  For  a  typical  two-dimensional  problem  having 
128  x  32  cells,  this  reduction  amounts  to  approximately  50  percent  and  30  per¬ 
cent  for  memory  and  arithmetic  operations,  respectively.  An  added  benefit  of 
the  nr.  >dification  is  that  it  simplifies  the  task  of  coding  the  scheme. 

The  multi-grid  feature  described  and  employed  by  Nayani6  can  be  used  with 
this  scheme  thereby  speeding  convergence  in  the  same  manner  as  he  demon¬ 
strated  with  the  original  scheme. 
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APPENDIXES 


APPENDIX  1 


DERIVATION  OF  Cjyk 

The  derivation  of  Cjik  is  accomplished  heuristically  by  comparison  of  the 
two-dimensional  case  with  the  one-dimensional  case.  In  this  section,  that  com¬ 
parison  will  be  presented  to  justify  the  terms  in  Chk.  Note  that  Cjik  simply 
multiplies  dQjtk.  Therefore,  approximations  in  this  coefficient  will  not  affect  the 
accuracy  of  the  results  since  the  terms  of  Resjtk  are  never  tampered  with. 

In  deriving  Equation  (36),  Nayani6  has  shown  an  earlier  st?p  to.be 


C  LjtkAQ  j— i>jt  +  CRjlkAQj+ltk  +  CC jikAQ  j,k 


+  CBj,k&Qj,k-\  +  CTj>k^Qj,k+i 

7-1  7-1 

_  D,  n  J  —  j  ’k  r(4)  'k  r(4)  A/^ 

-  &esj,k  a r  j-^kAQj~2'k  ~  At  }+±,kAQj+2'* 


7-1  7-1 

-(4)  Ar)  .L  _J± tl,(4)  , 

~~A7~  j’k-k  A(3l,*-2-  Ar  eM+iA  Qhk+2 


(48) 


For  Uin  —  1,  Equation  (35)  becomes 


AQ(m+i)  =  A  Q(m)  +  dQ 

Consider  a  sweep  that  is  left-to- right  and  bottom- to- top  (i.e.,  j  increasing  and 
k  increasing). 


CCjikdQjtk  +  CRj,kdQj+i'k  +  CTj<kdQjik+i 

J-1 


r- 1 


Jj+hk 


+  =  MS  (49) 


(4) 


RHS  =  Res]  k  -  Lj<k(AQ) 


where  i?es"fc  is  shown  in  Equation  (37)  and  Ljyk{AQ)  is 

MAO)  =  CLj.tAQ^  +  +  CCjftAQ% J> 

+  CRj.tAQ^l,,  +  CTj.kAQ($+l 

J~l  7-1 

J-j.fc  (4)  An(m+1)  j+$,k  (4)  AO(m) 

+  Ar  Vi*  V>~2-*  +  Ar  7+i.*A(*W>* 

J-1  J-l 

i  J’*~2  r(4)  4-  r(4)  AD^m^ 

+  Ar  i,k-\^j,k-2  +  Ar  €j,fc+iACA,*+ 2 


Equation  (49)  is  analogous  to  the  one-dimensional  Equation  (9) 


duj  -(-  -(ditj+i)  +  r(duj+2  -  4duJ+1  +  6 duj)  =  RHS  (9) 

RHS  =  Res]  -  Lj(Au) 

Lj(Au)  =  Auf>  +  §(Au<™>  -  Au^r’) 

+  r( Au£>  -  4Au$+j  +  6AuJm)  -  4Au£+1)  +  Au£+1}) 

During  the  process  of  sweeping  through  the  flow  field,  the  left-hand  side  terms 
of  Equation  (49)  represent  unknowns,  while  the  right-hand  side  terms  represent 
known  quantities. 

The  first  approximation  for  the  one-dimensional  case  is 


duj+ 2  ss  duj+i  fa  duj 

Similarly,  the  first  approximation  for  the  two-dimensional  case  is 

dQ j+2,k  ~  dQ j+\tk  fa  dQj<k+2  ^  dQj<k+\  %  dQjtk 
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In  one  dimension,  Equation  (9)  becomes 


(1  +  |  +  3  r)duj  =  RHS 


(50) 


In  two  dimensions,  Equation  (49)  becomes 


(coef)dQ  jtk  =  RHS 


(51) 


where 


coef  =  CCjzk  +  CRj'k  +  CTjtk  + 


j-i 

+  (4) 

At  €j+bk 


1  + 


Ji*+h  rW  T 
At  i>*+2 


The  second  modification  for  the  inner  iteration  in  one  dimension  is  to  replace 
c  by  |c|  and  for  all  sweeps  maintain  the  coefficient  of  duj  as 


(1  +  Y  +  3r) 


The  corresponding  modification  in  two  dimensions  is  a  series  of  approximations. 
Expanding  the  coefficient  of  dQjtk  in  Equation  (51)  into  all  its  terms  gives 


j-i  j-i 

coef  =  CCj,k  +  CR,,k  +  ctm  + 

—  ( inviscid  terms )  +  ( viscous  terms)  +  ( dissipation  terms ) 


where 


(inviscid  terms)  =  (\  +  \)(Vr,An  -  xvBn)j+^k  -  \(yvAn  - 

+  (2  +  2^~V(An  +  X(Bn^j’k+2  ~  \(~y*An  +  xtBn)j,k-± 


J7k 

+  -ir-1 

At 


(viscous  terms)  —  Re  J[(l  —  l)Rj  fc+i  + 


J-1  J_1 

( dissipation  terms)  ==  (1  —  l)- H — "  ej-^ 

+  (1  -  i)i±tie(2)  /  +  i±dc(.2)  x/ 

v  7  At  J>fc+j  At  ’ fc  - 


~j,k- 


+  (1  +  3  —  3) 
+  (1  +  3  -  3) 


7_1 

i+£,*  (4) 

At  i+i>* 

i±tie(4) 

At  i-fc+j 


/  +  (3  —  1) 
I  +  (3-l). 


(4)  f 

At  i-bk 

Jj,k-$  (4)  f 

Ar  i.*-i 


The  inviscid  terms  at  the  half  cell  locations  are  approximately  the  same  as  at 
the  cell  centers.  That  is 


(yvAn  -  xvBn)j+^k  ~  (yvAn  -  xnBn)j,k 
( yvAn  -  xvBn)j_L  k  ~  ( yvAn  -  xvBn)itk 
(-y<An  +  X(Bn)jMi  ~  (~y(An  +  x(Bn)jtk 


(-yfAn  +  x(Bn)j  k_x  ~  (~y^An  +  x^Bn)hk 


This  gives 


(inviscid  terms)  =  -(y^A™ 
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The  spectral  radius  of  the  first  two  terms  of  the  the  inviscid  terms  is  taken  next. 


SpRad[^(yvAn  -  Xr,Bn)j<k  +  ^(-y€An  +  x?Bn)j)fc] 

SpRad(y,jAn  -  xnBn)jtk  +  ^SpRad(-y^An  +  x(Bn)jik 

-\[\yr,u  -  x,u|  +  (x$  +  y2)*a]  +  i[|  -  y^u  +  x€v|  +  (x\  4-  yf)*a] 

where  u  and  v  are  velocities  and  a  is  the  speed  of  sound.  CFL  for  the  two- 
dimensional  case  is  defined  as 

_  ( velocity)(Atime )  [( velodty)(Adistance)](Atime ) 

(A  distance)  [area) 


where 


[(ue/ocity)(  Advance)]  =  [|y,u  -  xvv\  +  (x*  +  y*)^a]  +  [|  -  y5u  +  x^u|  +  (x|  +  y|)*a] 
A  time  =  At 
area  =  J~£ 

At  each  grid  node,  a  local  time  step  is  estimated  which  keeps  the  local  Courant 
number  approximately  constant.  That  estimate  is  given  by 


At  = 


[\llvu 


_ (CFL){J~l) _ 

xvv\  +  (xjj  +  y*)$a]  +  [|  -  y£u  +  x{u|  +  (x|  +  y| ) 2 <2 
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The  approximation  to  the  inviscid  terms  with  local  time-stepping  can  be  rewrit¬ 
ten  as 


1  J~l 

(inviscid  terms )  ~  -CFL  ?'  I  + 

2  A  r 


where  Ar  is  the  local  time  step. 

To  make  the  viscous  terms  insensitive  to  sweep  direction,  the  following  ap¬ 
proximation  is  made. 

( viscous  terms)  =  i?e-1[(l  —  l)Rj  ,*+*  + 

^^Re-\RjikH+Rjtk_^) 

Finally,  to  make  the  dissipation  terms  insensitive  to  sweep  direction,  the 
following  approximations  are  made. 


J  —  1  J  —  1  7  —  1  7—1 

(1  11  j  |  (2)  r  t  |  i~h'k  (2)  rs. 

{  Ar  i+hkI+  Ar  “  2^  Ar  >+*+/  +  Ar 

J-’ 


(1-1)- 


2 zidJ2) 

Ai 

-l 


(1+3-3) 

7: 

(l  +  3_3)-2£± 
Ar 


r— 1  7— 1  t — 1  t— l 

J'k+j  e(2)  ,/+-^ie(2)  ,  J~  If^±tie(2>  J  |  *fc-*c(2>  n 

Ar  i-*+s  Ar  >>fc~s  2^  Ar  >-fc+i  Ar  >.*-i 

■*—  1  T  — 1  T  — 1 

J  i-X 
J  2  . 


Ji±i+e(4)i  /  +  3(Vt+e(4)  7,  izL£e<«> 

Ar  '  Ar  >_£-*  2^  Ar  >+£>*  Ar  + 

r-i  7-1  j-l  7-1 

i±tLc(4)  I  I  11  J,*~*r(4)  3- f  _>»*+*  g(4)  7  |  >■*-*-( 4) 

Ar  >>*+i  Ar  2  Ar  >++i  Ar  i.*“i 


/) 

0 


The  resulting  approximation  for  the  dissipation  terms  then  becomes 


7—1  7—1 

~  1  /_>±i+  ,(2)  r  |  >-j.*  (2) 

—  2  Ar  >+*•*  Ar  Vi* 
J-1  J-1 

,  j,k+  j  (2)  r  ,  J.*~2  (2)  r\ 

Ar  i«*+i  Ar  M— i 

I  Z,Jj+hk  (4)  j-  .  (4)  r 

+  2(  Ar  >+i*i  +  Ar  V*,*7 

+  i|±it(*)  r+fiLitw  /) 

Ar  J-*+T  Ar  **-* 


(> dissipation  terms )  ~  ^(-r^i~e}+^ 


er,  ./ 
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Combining  the  approximations  for  the  inviscid  terms,  the  viscous  terms, 
and  the  dissipation  terms  gives  Cjtk- 


C„t  =  CFL(%i)( I  +  A-)l  +  *[%!«&  J 


I'1 


+  zJi±i'k e[4) i  j+ 


2  Ar  J+i.* 


Ar  T  Ar  £i-hl 

J-1  J-1 

J.*+2  (2)  .  o  J.fc+ 2  (4)  j- 

+  At  At  M+K 

i/  +  3^^ie(4.)  ,/ 
2 


1  +  3- 


Af  >-*.* 


+ 


J-,1  , 


Ar 


Ar 


(38) 
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APPENDIX  2 


* 


ADDITIONAL  SMOOTHING  FACTORS  IN  SCHEME 

To  make  the  original  scheme  represented  by  Equations  (34)-(36)  more  ro¬ 
bust,  additional  smoothing  factors  have  been  added  to  the  iteration  part  of  the 
scheme.  They  appear  in  Nayani’s  original  code,  but  have  not  been  reported  on 
in  his  dissertation6.  These  smoothing  factors  are  used  in  the  modified  scheme 
as  well.  The  purpose  of  this  section  is  to  derive  these  factors.  Note  that  these 
factors  affect  only  the  AQ  or  dQ  terms  and  in  no  way  alter  the  Res™ k  term,  so 
accuracy  of  the  scheme  is  not  compromised  by  this  addition. 

The  thin-layer  Navier-Stokes  equations  in  nondimensional  form  for  general¬ 
ized  curvilinear  coordinates  was  given  in  Section  3.1  as 


IrtJ-’O)  +  ~(,y,E  -  x,F) 


+  °-{-y(E  +  xiF)-Re-As  =  0 


(33) 


After  spatially  integrating  Equation  (33)  using  a  finite-volume  approach,  an 
Euler  implicit  scheme  is  used.  The  delta  form  of  the  equation  becomes 


-£;C>Qlt  +  Kv„A’  -  *,B")A<2"]  ^  „  +  (( -y(A"  +  *fB")A<?"] 


i-±,k 


i.*- 2 


-Re~\ASn) 


-  Ad 

-  Ad 

3k+*  =  #es"fc  (52) 


where 


dQ’  dQ 


and  Res™ k  is  given  by  Equation  (37). 
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Consider  the  second  term  of  Equation  (52).  It  is  a  central  difference  repre¬ 
sentation  for  the  derivative  of  the  bracketed  term. 

[( V,A"  -  x,B")AQ"\  ~  -  x„B»)A<?"]  (53) 

]-bk 

The  derivative  in  Equation  (53)  can  be  split  into  two  parts  (dropping  the  time 
level  n  for  convenience) 

A  A  A 

~[(y,A  -  x„B)AQ]  =  -g-^AAQ)  -  ^[x„BAQ]  (54) 

The  first  term  on  the  right  side  of  Equation  (54)  will  serve  as  a  pattern  for 
the  other  derivatives.  Instead  of  using  central  difference  to  represent 
an  upwind  scheme  is  used.  For  an  upwind  scheme,  A  is  expressed  as  follows.  ( B 
would  be  expressed  similarly.) 

A  =  A+  +  A~ 

A+  =  A  +  \A\  A-  =  A~\A\ 

2  ’  2 

where  |x4|  is  the  spectral  radius  of  A  times  the  identity  matrix. 

|A|  =  SpRad(A)I ,  SpRad(A)  =  |tt|  +  a 

u  is  velocity  and  a  is  the  speed  of  sound.  has  non-negative  eigenvalues  and 
A~  has  non-positive  eigenvalues. 

The  upwind  scheme  uses  backward  difference  for  A+  and  forward  difference 
for  A~.  The  difference  form  of  the  derivative  in  Equation  (54)  becomes 

A 

^(y,AAQ)  =  +  4-)AQ] 


(55) 
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J 


4 


J^(!/,AA<3)  ~  (v,A+AO)m  -  (M+A(3)Hl, 

+  (y,A-AQ)J+i,it  -  (yr,A~ AQ)jtk  (56) 

Substituting  Equation  (55)  into  Equation  (56)  gives 

^(v,AAQ)  ~  i[(v,(X  +  |A|)A(J)j,t  -  (»,(A  +  lADAO^-,,*] 

+  ~  |A|)AQ)>+i,t  -  (»,(A  -  W)AQ)Jl4] 

~  |[(yIJAAQ)i+1)*  -  (j/^AQ)^!,*] 

-  i[(»,|^|A(J)i+,,»  -  2(j,,|^|AQ)m  +  (j/ijI^|AQ)j_Ii4]  (57) 

The  first  bracketed  term  of  Equation  (57)  is  simply  the  central  difference 
representation  of  ^(y,«AAQ),  that  is 

\[(yr,AAQ)j+hk  -  (yjlAAQ)j-itk\  =  (y,  AAG)|£|’J 

The  second  bracketed  term  of  Equation  (57)  is  the  additional  smoothing  that 
results  from  using  upwind  differencing  on  the  AQ  side  of  Equation  (52).  It  is 
the  central  difference  representation  for  the  second  derivative  for  (y^AAQ),  that 
is 

-5fe,WA<?W  -  2(»,|X|A Q)j,t  +  (y,M|A <?)>_,.»]  =  -i«|(Sr,|A|A«)j,t 

where  6 1  is  the  central  difference  operator  representing  jp. 

The  pattern  is  the  same  for  -^(x^BAQ),  -^(—y^AAQ),  and  -^{x^BAQ). 
When  upwind  differencing  is  accomplished  for  all  of  these  terms,  then  Equation 
(52)  becomes 


43 


r  —  1 


-jrt^Qlk  +  KM"  -  x,B")AQ")  ' '  *"  +  [(-y{A"  +  X(Bn)AQn] 


At 


J-$,k 


j,k+k 


-  (f°c){62c  [(y,|A"|  -  *„|B”|)AQ"1M  +  6 \  [(-y{|A*|  +  x{|.B"|)Ae”]M} 


-  Be”1  (AS”) 


i»*+i 


Jjk  2 


-  Ad 


J+j  .* 


k 

J  2>k 


-  Ad 


M+4 


i.fc-t 


=  Soji 


(58) 


The  (/oc)  coefficient  before  the  second  differences  gives  control  of  the  amount 
of  smoothing  that  is  used.  A  typical  value  is  ( fac )  =  0.1.  Equation  (58)  reflects 
the  equation  that  Nayani’s  code6  actually  solves,  not  Equation  (52).  In  the 
present  modified  code  also,  implicit  smoothing  terms  are  used  in  the  iteration 
process  with  (fac)  =  0.1. 
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